#ifndef MYMATH_H
#define MYMATH_H
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef WIN32
#include <xmmintrin.h>
#include <emmintrin.h>
#endif
struct Vector3f
{
	float x, y, z;
	const Vector3f operator* (const float &other)  const
	{
		Vector3f r;
		r.x = this->x * other;
		r.y = this->y * other;
		r.z = this->z * other;
		return r;
	}
	const Vector3f operator*= (const float &other)
	{
		this->x = this->x * other;
		this->y = this->y * other;
		this->z = this->z * other;
		return *this;
	}
};
struct Vector4f
{
	float x, y, z, w;
};
class Matrix
{
public:
#ifdef WIN32
	__declspec(align(16)) float data[4][4];
#else
	float data[4][4];
#endif
	//float *data = new float[16];
	//data[0][0] data[1][0] data[2][0] data[3][0]
	//data[0][1] data[1][1] data[2][1] data[3][1]
	//data[0][2] data[1][2] data[2][2] data[3][2]
	//data[0][3] data[1][3] data[2][3] data[3][3]
public:
	operator float*()
	{
		return data[0];
	}
	void SetIdentity()
	{
		memset(data,0,sizeof(data));
		data[0][0] = data[1][1] = data[2][2] = data[3][3] = 1.0;
	}
	void SetOrtho(float left, float right, float bottom, float top, float zNear, float zFar)
	{
		memset(data,0,sizeof(data));
		data[0][0] = 2.0/(right - left);
		data[1][1] = 2.0/(top - bottom);
		data[2][2] = -2.0/(zFar - zNear);
		data[3][3] = 1.0;
		data[3][0] = -(right + left)/(right - left);
		data[3][1] = -(top + bottom)/(top - bottom);
		data[3][2] = -(zFar + zNear)/(zFar - zNear);
	}
	void SetPerspective(float fovyDegree, float ratio, float zNear, float zFar)
	{
		float f = 1.0/tan(fovyDegree/2);
		memset(data,0,sizeof(data));
		data[0][0] = f/ratio;
		data[1][1] = f;
		data[2][2] = (zNear + zFar)/(zNear - zFar);
		data[2][3] = -1.0;
		data[3][2] = (2 * zFar * zNear)/(zNear - zFar);
	}
	void SetNormalMatrix(const Matrix * modelView)
	{
		//NormalMatrix is 3x3 matrix
	}
	void SetTranslate(float x, float y, float z)
	{
		SetIdentity();
		data[3][0] = x;
		data[3][1] = y;
		data[3][2] = z;
	}
	void SetRotateX(float radian)
	{
		SetIdentity();
		data[1][1] = cos(radian);
		data[2][2] = cos(radian);
		data[1][2] = sin(radian);
		data[2][1] = -sin(radian);
	}
	void SetRotateY(float radian)
	{
		SetIdentity();
		data[0][0] = cos(radian);
		data[2][0] = sin(radian);
		data[0][2] = -sin(radian);
		data[2][2] = cos(radian);
	}
	void SetScaleUniform(float scale)
	{
		data[0][0] = scale;
		data[1][1] = scale;
		data[2][2] = scale;
		data[3][3] = 1.0f;
	}
	/*
	void mmul_sse_dot(const float * a, const float * b, float * r)
	{
	  __m128 a_line, b_line, r_line;
	  float mc[16] __attribute__((aligned(16)));  // 16-byte aligned temp array
	  for (int i=0; i<16; i+=4) {
		b_line = _mm_load_ps(&b[i]);              // b_line = vec4(column(b, i))
												  // remember that we are column-major
		for (int j=0; j<4; j++) {
		  a_line = _mm_set_ps(a[j+12], a[j+8], a[j+4], a[j]); 
												  // a_line = vec4(row(a, j))
												  // note that SSE is little-endian
		  r_line = _mm_mul_ps(a_line, b_line);    // r_line = a_line * b_line
		  _mm_store_ps(mc, r_line);               // copy r_line to memory
		  r[i+j] = mc[0]+mc[1]+mc[2]+mc[3];       // r[i][j] = sum(r_line)
												  //         = dot(a_line, b_line)
		}
	  }
	}
	*/
	const Matrix operator* (const Matrix &other)  const
	{
		Matrix r;
		for(int i = 0;i<4;i++)
		{
			for(int j = 0;j<4;j++)
			{
				r.data[i][j] = this->data[0][j]*other.data[i][0] + this->data[1][j]*other.data[i][1] + this->data[2][j]*other.data[i][2] + this->data[3][j]*other.data[i][3];
			}
		}
		return r;
	}
	const Vector4f operator* (const Vector4f &other)  const
	{
		Vector4f r;
		r.x = this->data[0][0]*other.x + this->data[1][0]*other.y + this->data[2][0]*other.z + this->data[3][0]*other.w;
		r.y = this->data[0][1]*other.x + this->data[1][1]*other.y + this->data[2][1]*other.z + this->data[3][1]*other.w;
		r.z = this->data[0][2]*other.x + this->data[1][2]*other.y + this->data[2][2]*other.z + this->data[3][2]*other.w;
		r.w = this->data[0][3]*other.x + this->data[1][3]*other.y + this->data[2][3]*other.z + this->data[3][3]*other.w;

		return r;
	}
	static void MatrixMultiply(const Matrix * a, const Matrix * b, Matrix * r)
	{
		for(int i = 0;i<4;i++)
		{
			for(int j = 0;j<4;j++)
			{
				r->data[i][j] = a->data[0][j]*b->data[i][0] + a->data[1][j]*b->data[i][1] + a->data[2][j]*b->data[i][2] + a->data[3][j]*b->data[i][3];
			}
		}
	}
	static void MatrixMultiplySSE(const Matrix * a, const Matrix * b, Matrix * r)
	{

	}
};
Vector3f Normalize(const Vector3f &v);
Vector3f Cross(const Vector3f &u, const Vector3f &v);
float Dot(const Vector3f &u, const Vector3f &v);

#endif
/*
http://fhtr.blogspot.com/2010/02/4x4-float-matrix-multiplication-using.html
http://unspecified.wordpress.com/2012/06/21/calculating-the-gluperspective-matrix-and-other-opengl-matrix-maths/
Normal matrix http://www.lighthouse3d.com/tutorials/glsl-tutorial/the-normal-matrix/
Generate mipmaps http://www.g-truc.net/post-0256.html
http://www.geeks3d.com/20110405/fxaa-fast-approximate-anti-aliasing-demo-glsl-opengl-test-radeon-geforce/
https://developer.nvidia.com/csaa-coverage-sampling-antialiasing
*/